	* Estimate parametric CZ F.E. models on the same sample as the TWFE model

	use "${clean_data}/panel_physicians_clean.dta", clear	
	keep id year cz90 year logptotinc female age is_foreign race married pr_taxonomy_code birthplace_char ein
	fmerge 1:1 id cz90 year using "${output}/twfe/geo_twfe/physicians/IndividFixedEffects_ALL_physicians", keepusing(id) keep(match) nogen
	egen prtaxm_num=group(pr_taxonomy_code)

	* Estimation
	foreach x in RAW RESIDUALIZED {
		
		if "`x'" == "RAW" {
			local fes "cz90" 
			local labelfe "CZ"
		}
		else if "`x'" == "RESIDUALIZED" {
			local fes "cz90 year female age is_foreign race married prtaxm_num birthplace_char ein" 
			local labelfe "Residualized: Yr/F/A/US/R/M/PrTax/PO/EIN" 
		}

		di "Regress loginc on `labelfe'"
		su logptotinc
		reghdfe logptotinc, a(`fes', savefe) 
		
		cap gen sample = e(sample)
		
		loc i = 0
		foreach fe of loc fes {
			loc i = `i' + 1
			ren __hdfe`i'__ fe_`fe'_`x'
			label variable fe_`fe'_`x' "FEs resulting from regression of log(income) on `labelfe'-FEs"
		}
	}
	
	keep id year cz90 fe_cz90_RESIDUALIZED fe_cz90_RAW
	isid id year
	save ${temp}/parametric_CZFE.dta, replace 

	* Correlate raw, parametrically adjusted, and non-parametric CZ F.E. and individual FE
		
	use ${temp}/parametric_CZFE.dta, clear
	fmerge 1:1 id year cz90 using "${output}/twfe/geo_twfe/physicians/IndividFixedEffects_ALL_physicians", assert(match) /// 
		keepusing(__hdfePLACE__ __hdfeINDIVID__ logptotinc) nogen
		
	keep id cz90 __hdfeINDIVID__  __hdfePLACE__ fe_cz90_RESIDUALIZED fe_cz90_RAW
	duplicates drop
	egen aux=max(fe_cz90_RESIDUALIZED), by(cz90)
	replace fe_cz90_RESIDUALIZED=aux
	drop aux
	duplicates drop	


	foreach var in fe_cz90_RAW fe_cz90_RESIDUALIZED __hdfePLACE__ {		
		
		reg `var' __hdfeINDIVID__, vce(cluster cz90)
		gen rsquared=e(r2)
		matrix b=e(b)
		gen slope=b[1,1]
		gen cons = b[1,2]
		matrix se=e(V)
		gen se=sqrt(se[1,1])
		binscatter  `var' __hdfeINDIVID__, xtitle("Individual FE from TWFE model", height(10)) ytitle("`var'") genxq(bin)
		
		gegen idfe=mean(__hdfeINDIVID__), by(bin)
		gegen czfe=mean(`var'), by(bin)
		gegen n_ur=count(id), by(bin)
		
		preserve
		keep bin idfe czfe slope cons se rsquared n_ur 
		duplicates drop
		drop if missing(bin)
		sort bin
		gen specification="`var'"
		save ${temp}/`var'.dta, replace
		restore 
		
		drop bin idfe czfe slope cons se rsquared n_ur
	}  
	
	* Create a version of binscatters that follows Card et al approach and collapses
	* person F.E. to CZ level before measuring the correlation
	
	use ${temp}/parametric_CZFE.dta, clear
	fmerge 1:1 id year cz90 using "${output}/twfe/geo_twfe/physicians/IndividFixedEffects_ALL_physicians", assert(match) /// 
	keepusing(__hdfePLACE__ __hdfeINDIVID__ logptotinc) nogen
		
	keep id cz90 __hdfeINDIVID__ __hdfePLACE__ fe_cz90_RESIDUALIZED fe_cz90_RAW
	duplicates drop
	egen aux=max(fe_cz90_RESIDUALIZED), by(cz90)
	replace fe_cz90_RESIDUALIZED=aux
	drop aux
	duplicates drop
	egen mean__hdfeINDIVID__=mean(__hdfeINDIVID__), by(cz90)
	keep cz90 mean__hdfeINDIVID__ __hdfePLACE__ fe_cz90_RESIDUALIZED fe_cz90_RAW
	duplicates drop
	isid cz90
		
	foreach var in fe_cz90_RAW fe_cz90_RESIDUALIZED __hdfePLACE__ {		
		
		reg `var' mean__hdfeINDIVID__, robust
		gen rsquared=e(r2)
		matrix b=e(b)
		gen slope=b[1,1]
		gen cons = b[1,2]
		matrix se=e(V)
		gen se=sqrt(se[1,1])
		binscatter  `var' mean__hdfeINDIVID__, xtitle("Mean individual FE from TWFE model", height(10)) ytitle("`var'") genxq(bin)
		
		gegen idfe=mean(mean__hdfeINDIVID__), by(bin)
		gegen czfe=mean(`var'), by(bin)
		gegen n_ur=count(cz90), by(bin)
		
		preserve
		keep bin idfe czfe slope cons se rsquared n_ur 
		duplicates drop
		drop if missing(bin)
		sort bin
		gen specification="`var'"+" ID FE averaged"
		save ${temp}/`var'_idcollapsed.dta, replace
		restore 
		
		drop bin idfe czfe slope cons se rsquared n_ur
	}  
	
	use ${temp}/fe_cz90_RAW.dta, clear	
	append using ${temp}/fe_cz90_RESIDUALIZED.dta
	append using ${temp}/__hdfePLACE__.dta
	append using ${temp}/fe_cz90_RAW_idcollapsed.dta
	append using ${temp}/fe_cz90_RESIDUALIZED_idcollapsed.dta
	append using ${temp}/__hdfePLACE___idcollapsed.dta
	
	order specification bin czfe idfe slope se rsquared n_ur
	isid  specification bin
	
	drbest_docinc czfe idfe slope se rsquared  
	drbcount n_ur, gen(n)
	order specification bin czfe idfe slope se rsquared n_ur n

	export delimited using "${mypath}/intermediate_csv/geotwfe0405-indivfe_vs_czfe_manualczfe.csv", replace dataf delim(tab) 	
